#Replication package for Kim, DellaPosta, and Essig
#Code for reproducing agent-based model

rm(list=ls())

library(purrr)
library(dplyr)
library(MASS)
library(clusterGeneration)

#######################################################
#FUNCTIONS TO CALL LATER 

impute_positions = function(PolID, VotID, PublicIssues, Politicians, VoterAssoc, AllIssues, Continuous, Slope){  
  
  #Weight known positions by comparing with party
  IssueWeights = sapply(PublicIssues, function(iss) {
    weight = (Politicians[PolID, iss] - mean(VoterAssoc[[VotID]][,iss], na.rm=T))/2
    return(as.numeric(weight))
  })
  
  #Impute unknown positions
  IdeolScore = mean(IssueWeights)
  Positions = sapply(1:AllIssues, function(iss) {
    if(iss %in% PublicIssues){
      p = Politicians[PolID, iss]
    }else{
      imp_pos  = ifelse(!(is.na(mean(VoterAssoc[[VotID]][,iss], na.rm=T))),
                     (mean(VoterAssoc[[VotID]][,iss], na.rm=T)+IdeolScore),
                     IdeolScore)
      if(Continuous==T){
        p = imp_pos
        if(p>1){
          p=1
        }
        if(p<-1){
          p=-1
        }
      }else{
        posprob = 1/(1+exp(-1*Slope*(imp_pos)))
        p = ifelse(rbernoulli(1, p=posprob)==T, 1, -1)
      }
    }
    return(p)
  })
  return(Positions)
}

#########################################3
#CODE FOR SIMULATION

party_model = function(n_voters=10000, #voter agents
                       n_pols=50, #politician agents
                       k=20, #issue dimensions
                       h=3, #visible positions
                       continuous_positions=F, #for "nuanced positions" condition, change to true
                       unequal_weight=F, #for "varying salience" condition, change to true
                       biased_display=F, #for "biased display" position, change to true
                       slope=25, #beta parameter
                       measurement_interval=1000, 
                       stop_cor=.99, #how strong between-politician correlations must be for model to stop
                       agent_memory=10, #m parameter 
                       voter_partisanship = 2){ #x parameter
  
  #Create politician population
  if(continuous_positions==T){
    pols = data.frame(matrix(runif(n_pols*k, min=-1, max=1), nrow = n_pols, ncol = k))
  }else{
    pols = data.frame(matrix(sample(c(-1, 1), n_pols*k, replace=T), nrow = n_pols, ncol = k))
  }
  pols$id = 1:n_pols

  #Create voter population
  issue_strengths = runif(k, min = -1, max = voter_partisanship)
  issue_corrs = rcorrmatrix(k, alphad = 1)
  voters = mvrnorm(n = n_voters, issue_strengths, issue_corrs)
  voters = as.data.frame(voters)

  #voter initially randomly supports pol from their party
  voters$pol = sample(pols$id, n_voters, replace=T)
  
  #Initialize cognitive associations
  a = list()
  mat = matrix(NA, nrow = agent_memory, ncol = k) #each agent can remember a fixed number of observations
  for(y in 1:n_voters) {
    a[[y]] = mat
    }
  
  #MAIN BODY OF MODEL
  polcor = 0
  i = 0
  
  while(polcor<stop_cor){
    i = i+1

    #Sample voter to update
    agent = sample(1:n_voters, 1)
    
    #Sample two politicians--one the voter currently supports, and a randomly chosen other
    CurrentPol = voters$pol[agent] #currently supported politician
    SecondPol = sample(pols$id[-CurrentPol], 1)
  
    #Politicians choose salient issues to display
    if(biased_display==T){
      if(h<k){
        known_issues = sample(1:k, h, replace=F,
                              prob = abs(colMeans(pols[,1:k]))/sum(abs(colMeans(pols[,1:k]))))
      }else{
        known_issues = 1:k
      }
   }else{
      known_issues = sample(1:k, h, replace=F)
   }
    
    #Voter updates associations
    for(m in known_issues){
      a[[agent]][,m]=lag(a[[agent]][,m], n=2)
      a[[agent]][1,m]=pols[CurrentPol,m]
      a[[agent]][2,m]=pols[SecondPol,m]
    }

    #Voter imputes pols' positions
    CurrentPolPos = impute_positions(CurrentPol, agent, known_issues, pols, a, k, continuous_positions, slope)
    SecondPolPos = impute_positions(SecondPol, agent, known_issues, pols, a, k, continuous_positions, slope)
  
    #Voter compares positions to their own
    if(unequal_weight==T){
      weights = abs(as.numeric(voters[agent,1:k]))/sum(abs(as.numeric(voters[agent,1:k])))
      CurrentPolDiff = sum(abs(as.numeric(voters[agent,1:k]) - CurrentPolPos)*weights)/k
      SecondPolDiff = sum(abs(as.numeric(voters[agent,1:k]) - SecondPolPos)*weights)/k
    }else{
      CurrentPolDiff = sum(abs(as.numeric(voters[agent,1:k]) - CurrentPolPos))/k
      SecondPolDiff = sum(abs(as.numeric(voters[agent,1:k]) - SecondPolPos))/k
    }

    #Voter chooses which pol to support as probabilistic function of match with each
    #Losing pol copies a position from winning pol
    CurrentPolProb=1/(1+exp(-1*slope*(SecondPolDiff-CurrentPolDiff)))
    updated_pos = sample(1:k, 1)
    if(rbernoulli(1, p = CurrentPolProb)==T){
      pols[SecondPol, updated_pos] = pols[CurrentPol, updated_pos]
      }else{
        voters$pol[agent]=SecondPol
        pols[CurrentPol, updated_pos] = pols[SecondPol, updated_pos]
        }
  
    #At fixed intervals, check for convergence on a single position set 
    if( (i/measurement_interval) == round(i/measurement_interval) | i == 1) {
    
      pos_cor_mat = matrix(NA, nrow=n_pols, ncol=n_pols)
      for(x in 1:n_pols) {
        for(y in 1:n_pols) {
          if(sum(as.numeric(pols[x,1:k])-as.numeric(pols[y,1:k]))==0){
            pos_cor_mat[x,y]=1}else{
              pos_cor_mat[x, y]=cor(as.numeric(pols[x,1:k]), as.numeric(pols[y,1:k]))
            }
          }
      }
      }
    pos_cor_mat[lower.tri(pos_cor_mat)] = NA #preference similarity only measured for unique pairs
    diag(pos_cor_mat)=NA #set diagonals to 0
    polcor = mean(pos_cor_mat, na.rm=T)
  }
  
  #record proportion of issues on which voters and politicians have pro-party positions on average
  voters_1 = length(which(colMeans(voters[,1:k])>0))/k
  pols_1 = length(which(colMeans(pols[,1:k])>0))/k
  
  outcome = paste(voters_1, pols_1, sep=",")
  return(outcome)
}
################################################
#Test out model in individual runs

r = party_model(n_voters=10000, 
                n_pols=50, 
                continuous_positions=F,
                biased_display=F,
                unequal_weight=F,
                k=20, 
                h=5,
                slope=25, 
                measurement_interval=1000,
                stop_cor=.99,
                agent_memory=10,
                voter_partisanship=2)

r

################################################
#Test conditions used to generate Figure 2

replications = 500 #set number of replications per condition

#Parameters to manipulate
public_issues = seq(2, 20, 2)
partisanship = c(1, 2, 3)

params = expand.grid(public_issues, partisanship)
names(params) = c("Public_Issues", "Voter_Partisanship")
runs = as.data.frame(params[rep(seq_len(nrow(params)), each = replications), ])

results = lapply(1:nrow(runs), function(t) {

  r = party_model(n_voters=10000, 
                  n_pols=50, 
                  continuous_positions=F,
                  biased_display=F,
                  unequal_weight=F,
                  k=20, 
                  h=runs$Public_Issues[t],
                  slope=25, 
                  measurement_interval=1000,
                  stop_cor=.99,
                  agent_memory=10,
                  voter_partisanship=runs$Voter_Partisanship[t])
  print(t)
  return(r)
})
#############################################
#Test alternative conditions from Appendix results

#BIASED DISPLAY
biased = lapply(1:replications, function(t) {
  r = party_model(n_voters=10000, 
                  n_pols=50, 
                  continuous_positions=F,
                  biased_display=T,
                  unequal_weight=F,
                  k=20, 
                  h=5,
                  slope=25, 
                  measurement_interval=1000,
                  stop_cor=.99,
                  agent_memory=10,
                  voter_partisanship=3)
  print(t)
  return(r)
})

#VARYING SALIENCE
unequal = lapply(1:replications, function(t) {
  r = party_model(n_voters=10000, 
                  n_pols=50, 
                  continuous_positions=F,
                  biased_display=F,
                  unequal_weight=T,
                  k=20, 
                  h=5,
                  slope=25, 
                  measurement_interval=1000,
                  stop_cor=.99,
                  agent_memory=10,
                  voter_partisanship=3)
  print(t)
  return(r)
})

#NUANCED POSITIONS
contin = lapply(1:replications, function(t) {
  r = party_model(n_voters=10000, 
                  n_pols=50, 
                  continuous_positions=T,
                  biased_display=F,
                  unequal_weight=F,
                  k=20, 
                  h=5,
                  slope=25, 
                  measurement_interval=1000,
                  stop_cor=.99,
                  agent_memory=10,
                  voter_partisanship=3)
  print(t)
  return(r)
})

#MODERATELY HIGH NOISE
noise1 = lapply(1:replications, function(t) {
  r = party_model(n_voters=10000, 
                  n_pols=50, 
                  continuous_positions=F,
                  biased_display=F,
                  unequal_weight=F,
                  k=20, 
                  h=5,
                  slope=10, 
                  measurement_interval=1000,
                  stop_cor=.99,
                  agent_memory=10,
                  voter_partisanship=3)
  print(t)
  return(r)
})

#HIGH NOISE
noise2 = lapply(1:replications, function(t) {
  r = party_model(n_voters=10000, 
                  n_pols=50, 
                  continuous_positions=F,
                  biased_display=F,
                  unequal_weight=F,
                  k=20, 
                  h=5,
                  slope=5, 
                  measurement_interval=1000,
                  stop_cor=.99,
                  agent_memory=10,
                  voter_partisanship=3)
  print(t)
  return(r)
})
